Statistical mechanics of Monte Carlo sampling and the sign problem. 
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Monte Carlo sampling of any system may be analyzed in terms of an associated glass model - a variant of the 
Random Energy Model - with, whenever there is a sign problem, complex fields. This model has three types 
of phases (liquid, frozen and 'chaotic'), as is characteristic of glass models with complex parameters. Only the 
liquid one yields the correct answers for the original problem, and the task is to design the simulation to stay 
inside it. The statistical convergence of the sampling to the correct expectation values may be studied in these 
terms, yielding a general lower bound for the computer time as a function of the free energy difference between 
the true system, and a reference one. In this way, importance-sampling strategies may be optimized. 



The sign problem appears when we are asked to compute a 
partition function of the form: 
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and W = Wu + iWi is not real. In classical statistical me- 
chanics s denotes a configuration in space and W — /3E is 
proportional to the energy, while in the quantum case, s are in 
space-time and W is the quantum action (and we refer to this 
as 'quantum Monte Carlo'). The problem is that a stochas- 
tic sampling cannot be performed with the complete W, since 
then the probabilities are no longer positive and real. 

An imaginary part of W arises naturally in quantum me- 
chanics with real time, in field theories like QCD with 
baryons, and in some of the most important systems in con- 
densed matter. In many cases, it is introduced when W is real 
and of the form W — Wo ~ bJ^C^ one needs a form 
linear in the Cq,: 
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The Hubbard-Stratonovich transformation above introduces 
an imaginary term Wi = \/\b\Yl ^aCa in the repulsive case 
6 < 0. 

Complex fields: Although the discussion in this paper 
is general, it is instructive to concentrate on a (very large) 
class of systems for which the sign problem takes a partic- 
ularly simple form. They are characterized by having a real 
action W^i, plus an imaginary term Wi = ihjAI{s), where 
M{s) is an integer-valued function of the configuration, and 
appear in several contexts: 

• 0-terms in Euclidean field theories and solid state. M{s) 
is a topological number associated with the configuration. The 
case hj = 6 = m, when the partition function is real, is par- 
ticularly interesting: the term Wi may change an otherwise 
gapped problem into a gapless one (see e.g. ifTllSl fTOl ) 

• Fermion systems such as the Hubbard model 
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The number operators rii-f , may be decoupled with a 
Hubbard-Stratonovich transformation (either with either con- 



tinuous or spin f3\ variables, and the remaining determi- 
nants evaluated. When the model is repulsive, these determi- 
nants pick up a minus sign counting fermion world-line cross- 
ings |4|, and in some cases M{s) can be shown to be a topo- 
logical number related to the Hubbard Stratonovich field s |5 |. 
• An alternative strategy to uncouple the number operators 
introduced by DeForcrand and Batrouni |6J uses a modified 
uncoupling with spins, and then M(s) is the sum over space 
and time of the number of "up" spins. 

In all these cases the sign problem arises because we have 
to calculate a partition function 
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with F{AI) a real, extensive quantity that may be computed 
with standard Monte Carlo. This is a deceptively simple ex- 
pression, since the sum (j3]l contains detailed cancellations that 
require an enormous precision in the F{M). 

Sampling: In the presence of sign problem, the alternative 
is to use only the real part in the Monte Carlo rule, and evalu- 
ate the expectation of a function O as: 
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where (•)r denotes an average evaluated with the Boltzmann 
weight of the real part Wr, a quantity easily calculable with 
Monte Carlo sampling. Somewhat more generally, we may 
apply reweighting. 
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where the average is associated with the energy —A + 
Wft., involving the suitably chosen importance sampling func- 
tion A. In practice, these expectation values are generated by 
running a Monte Carlo program with energy —A + Wn for a 
time T- If the correlation time of the Monte Carlo dynamics 
is T, this is essentially like using R independent samplings: 
R^T /t , which we parametrize as R = e''^ . 
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We may restate the problem as follows: we are given R = 
e^^ configurations Sa = {si, sr}, chosen independently 
with probability 
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where the sum in the denominator runs over all possible values 
the configuration s' my take. We estimate averages (j6]l as: 
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which is encoded in the generating (partition) function: 



^ ^-A(s^)-iWj{sa)+jO{s^) 
a=l 
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Note that we use the same sampled configurations for both 
numerator and denominator. 

Glasses: The connection with glasses arises when we ask: 
what is the expected result of this procedure? In order to an- 
swer this, we need to compute the expectation value of ([8]l, 
or, equivalently, the expectation of the logarithm of Z: this is 
a quenched average, a consequence of the random variables 
appearing in both numerator and denominator of ([8]l. We cal- 
culate then: 



ln|2(j)| = 
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Real and imaginary parts of expectation values are obtained 
from: 
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We have chosen to express everything in terms of the average 
of the logarithm of the norm of the partition function, because 
its phase in general fluctuates from realization to realization 
(|7][8l. Equations (j7]i and (j9]) define a disordered system, a 
generalization of the Random Energy Model |9 1. The quantity 
ln|Z(j)| is self-averaging in the thermodynamic limit. The 
result contains, as we shall see, different phases, depending on 
the parameters of the original problem and on the number of 
samples e'^'^ : 7 is now a phase parameter. Of these phases of 
the glassy model, only one coincides with the original model, 
the others give - even as N 00 - wrong results for the 
sampling. Monte Carlo sampling in the presence of the sign 
problem may be understood as the art of being in the only 
'good' phase. 

For clarity of presentation, we shall now specialize to the 
situation described by complex fields ([3]) and Q, and then 
mention how to proceed generally. We collect a number 
of samples R, classify them with their value M, obtained 
with probability P{M). The result is a set of M{M) , with 
X]7\/A/'(M) = R = e^''. Expectation values are obtained 
from (cfr 
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The N are random variables, and we wish to calculate the 
expected result for ( [T2| . We recognize a form of the Random 
Energy Model (REM) with imaginary inverse temperature hi, 
and M playing the role of energy. The input from the original 
system is through the (real) free energy densities F{M), for 
real M, the quantity we estimate through our simulation. (The 
usual REM corresponds to a quadratic F{M)). The solution 
of this problem follows closely refs. fTHSj. 

For large N, because of independence, one can convince 
oneself that 

((AA(Af) - {U{M))f) ^ {Af{M)) ^ R P{M) (13) 

We have that either R P{M) <C 1, and there are no samples 
of that value of M, or R P{M) > 1. In the latter case, 
one has that the average number of samples of given A/ is 
M{M) ^ RP{M), but with an eiTor ^M{M), where: 
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For convenience we have defined the intensive quantities 
/(to = M/N) = F{M)/N. The term with /(to,,) comes 
from the normalization integral J dme^^^^^^K and TOq is the 
real axis saddle ^ (toq) = 0. 

In order to calculate the result of (12 1 we start with the par- 
tition function: 
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Writing Af - RP{M) + R^P^ {M)r]{M), where r]{M) is a 
random number of zero mean and unit variance, we may split 
Z into a deterministic Zdet plus a fluctuating part Z^rr- The 
deterministic part reads: 
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The integral (16 1 is dominated by the saddle point 
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Note that, unlike the real-axis saddle TOq, the total one rrisp is 
complex. This is the correct solution, and we shall denote it ' 
Phase r. 

The shifting of the correct saddle point into the complex 
plane is due to the rapid oscillations introduced by the imagi- 
nary term, which have the effect of making the final result no 
longer dominated by the envelope maximum at TOq - and in- 
deed, the final result is exponentially smaller than the peak 
g-Npf(ma} ^ In the case of our sampling, this means that 
we should also consider two effects, each one leading to a 
new phase: The first is the fact that for values of to dis- 
tant enough from the maximum, the sampling produces no 
configurations at all. The region of integration is restricted 
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to values m bounded by the points ni}, such that M{M) = 
g-/3Ar[/(m,)-/(m„)]^ > ^^lat is, within the Umits given by: 



as above, we conclude that the free energies differences now 
become: 



/3[/(mb) - /(mo)] = 7 



(18) 



The second effect comes from the fact that our integrand is 
noisy: this is particularly serious because the correct result 
depends upon very strict cancellations giving an exponentially 
smaller answer The error due to the noise of the integrand is: 



Zerr = J dm ^M{M)e-^'^''^ T]{m) 



This is a fluctuating quantity. The expectation of jZerr P 
over noise is easy to evaluate by saddle point, using 

{r]{m)rj{m')) ^ 5{m — m'), and one gets that: 
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All in all, there are three types of phases |7, 8|, a charater- 
istic feature of glassy systems with complex parameters: 

• Phase I: The correct sampling of the original system, (up 
to usual finite-size effects). Here it appears as a 'liquid' phase 
of the glass, dominated by a saddle point mgp that is outside 
the domain of integration of ( [T7| , either because m^p is com- 
plex, or because it is real but beyond the original interval. 

• Phase II: A 'frozen' phase, where the solution is domi- 
nated by one of the boundaries in sampling mb satisfying ( [T8] l. 

• Phase III: An 'chaotic' phase: the sampling noise destroys 
the cancellations brought about by the oscillatory terms, and 
the peak of the 'envelope' of the integrand dominates. This 
phase can be shown to contain a surface density of Lee- Yang 
zeroes in the complex plane \TJ^. 

In the present case, with A — Q, phase 11 never dominates. 
The boundary between Phases I and III is given by the line 
where the two real parts of the free energies coincide and 
when the contour deformation over the integral ( [T6] l passes 
through the saddle. Then, in order to be in the Phase I, we 
need 

hi max {mjy) > Im[/3/(msp) + ihimsp] > hj min (mb) 

(20) 

-[Re[/3/(m3p) + z/i/m,J - /Sfim,)] + ^ > (21) 

(note that f{msp) < f{mo)), which may also be written in 
terms of the free energy at hj and the free energy at hj = 0. 
This is the condition for the simulation to yield the good re- 
sult, which must be read as a bound on the simulation time 
through the condition that 7 be large enough. It incorporates 
the well-known exponential dependence on the inverse tem- 
perature, and a free energy difference between the true model 
and one without the sign problem (i.e. with real field). Note 
that increasing the system size N at fixed simulation time T 
is dangerous, as it implies making 7 smaller. 

Consider next introducing an importance-sampling func- 
tion of the form A = Na{m). Following the same reasoning 
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where to+ and m^ minimize /3/(m) ± a{m), respectively, 
and TOfc are given by the generalization of ([T8]i: 



7 - Pifimb) - fim-)) + a{mb) - a(m-) = (22) 

Boundaries and free energy differences have changed, and 
saddle may even be outside the domain of sampling between 
boundaries. All three phases are possible. 

The important point is that this phase structure, and hence 
the values of 7 for correct sampling, depends on the free en- 
ergy of the true system (which we cannot change), and on the 
free energies f{m^), f{mi,). The latter are properties of the 
system under real fields, and hence easy to calculate for any 
A. One may then choose the optimal function A that will min- 
imize the 7 needed to be in phase /. 

The general case: It is easy to extend the analysis of the 
disordered model for the general case, in which the imaginary 
part Wi takes an abitrary form. We have to express free en- 
ergies and probabilities in terms of three collective variables 
Wi and A, playing the role of M. The probability of a pair 
(Wi,A) is: 

PiWi.A) cx J2 e^'^'^-^^'-'^diWiis)-WiI))SiA{s)-A) 

s 

and we use it to compute the average of In \ Z{Wi, A)\. There 
may be again three phases, depending on whether what dom- 
inates is a saddle point of complex {Wi, A), the boundary in 
the (Wi, yl)-plane where there is sampling, or the incoherent 
noise. The frontiers of these phases depend on number of sam- 
ples through 7, and also on the choice of importance sampling 
function A. 

The Lee- Yang zeroes: At the Lee- Yang zeroes of the orig- 
inal problem (i..e. Phase T), Z = and / — > 00. It may seem 
that these are then problematic regions for sampling. Actually, 
this is not really so. The reason is as follows: in the large N 
limit, Lee-Yang zeroes appear when two solutions compete: 



-Nuh+ifnih) 



(23) 



(note that here we are using the Legendre-transformed free en- 
ergy, in terms of h). Along the transition line in the complex 
plane where f^{h) = f^{h) — fnih), the partition func- 
tion is Z ^ e-^/«(^)[e^^/^('') - e^^-^^C*)], and zeroes occur 
where the oscillating factors cancel. In order for the oscil- 
lating factor to be of the same order as Nfji{h), the factor 
[e^^ti W — e^^^' ^'''J has to be exponentially small in N, and 
this happens only very close to each zero. We conclude that 
each zero has a very small region of influence in which the 
partition function is really small. 
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FIG. 1 : Simulated minus exact value of the magnetization in terms 
of the complex fields /i for 7 — 0.5, J — 0.5 and A'^ = 100. The 
boundary between phases I and III is indicated with a dashed line 
determined through l |21| l. 



An example: ferromagnet in a complex field: Let us 

consider, as a simple example, a mean-field Ising ferromag- 
net under a complex field h ~ hu + ihj. We have that 

Z = Q-Nl3f(m)-ihimN ^ with: 
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(1 + m)ln(l + m) ~ -^C^ 
- hum 



m) ln(l — to) 
(24) 



tical mechanics of a glass model at complex fields. A first 
suggestion of this analysis is that simulations should be char- 
acterized by the size N and 7 — ^-j^, where T is the time. 
A first simulation with relatively modest values of N may be 
used to obtain the boundary jc of phase /, and only then em- 
bark on larger sizes, with times determined by fixed 7 ^ 7^. 

A more sophisticated strategy is to approach the desired pa- 
rameter (say, hi — tt) gradually, from smaller hj, and es- 
timate ^c{hi)- Because the sign problem is less severe for 
smaller values of hi, one may reach the desired situation by a 
succession of correctly sampled problems. 

As mentioned above, all importance sampling strategies 
work by increasing the free energy of phases II and III, thus 
making phase I more favorable at given value of 7. Because 
of the nature of these phases, this may be optimized on the 
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The correct solution of the problem is given by the saddle 
point equation to^^ = tanh(/3 Jw^p + h) 

For complex h, in general m^p will be complex, except for 
hi = ^, where it will be real but outside the interval [—1,1]. 

Figure [T] shows the sampled magnetization (subtracting, for 
clarity, the exact value) in terms of the field h, with a fixed 
sampling time e^''. In Figure|2]we show the frontier in detail, 
with the regions in Phase I surrounding the lee- Yang zeroes. 
Clearly the size of the islands decreases for larger systems. 

Conclusions: We have shown that the error brought about 
by the sign problem may be discussed in terms of the statis- 



FIG. 2: A detail of the frontier between phases I and III under the 
same conditions of Figure[T]with A'^ = 60 (left) and A'^ — 100 (right). 
The islands visible within phase I (correct sampling) are the neigh- 
borhoods of Lee- Yang zeroes of the ferromagnet 



basis of the knowledge of the free energy of the system at real 
fields alone. 

To conclude, let us mention that an analysis very similar to 
the one made here may be applied to the numerical implemen- 
tation of analytic continuation. 
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